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We consider the geometrical optics problem of finding a system of two reflectors that 
transform a spherical wavcfront into a beam of parallel rays with prescribed intensity 
distribution. Using techniques from optimal transportation theory, it has been shown 
before that this problem is equivalent to an infinite dimensional linear programming (LP) 
problem. We investigate techniques for constructing the two reflectors numerically. A 
straightforward discretization of this problem has the disadvantage that the number of 
constraints increases rapidly with the mesh size. So with this technique only very coarse 
meshes are practical. To address this well-known issue we propose an iterative solution 
scheme. In each step an LP problem is solved. Information from the previous iteration 
step is used to reduce the number of constraints necessary. As a proof of concept we 
apply our proposed scheme to solve a problem with synthetic data. We give evidence that 
the scheme converges. We also show that it allows for much finer meshes than a simple 
discretization scheme. There exists a growing literature for the application of optimal 
transportation theory to other beam shaping problems. Our proposed scheme is easy to 
adapt for these problems as well. 

1 Introduction 

The following beam shaping problem from geometric optics was described in [15]; see Figure 1: 
Suppose a point source emits a spherical wavefront with a given intensity distribution. The 
problem at hand consists of transforming this input beam into an output beam of parallel 
light rays with a prescribed intensity distribution. This transformation is to be achieved with 
a system of two reflectors. The problem has some practical importance in engineering, see 
further literature cited in [15]. 



The first rigorous mathematical solution to the problem was provided in [10], using an 
approach based on the theory of optimal transportation [5, 6, 8, 19]. See also the references 
[11, 9, 21, 12, 16] which deal with other beam shaping problems using related techniques. 

The result of [10] is summarized in section 2 below, with Theorem 2.5 stating the main 
result. The central feature is that the original reflector design problem is reformulated as an 
infinite dimensional constrained optimization problem, namely the problem of minimizing a 
certain functional on a function space. It is the dual problem for the problem of finding a map 
from the input aperture to the output aperture which minimizes a certain cost functional. (See 
Problem II in section 2 for the exact statement.) 
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This reformulation of the problem is not only of theoretical value for questions of existence 
and uniqueness of solutions, but it also translates into a practical method for finding the 
solution. In fact, the discretization of the infinite dimensional constrained optimization problem 
is a standard linear programming problem and can be solved numerically. In this paper, we 
describe a numerical scheme for solving these problems. An immediate obstacle is that the 
number of constraints becomes very large even for relatively coarse meshes if there arc M 
sample points in the input aperture (denoted by D in Figure 1) and points in the output 
aperture (T in Figure 1), then the number of constraints is M-N. We found that with a simple 
straightforward discretization method, standard LP solvers (on PCs with up to 4 GB RAM) 
could not handle more than approximately 500 sample points on the domains T and D - a 
number which is arguably too small for many applications. For this reason, we devised a more 
elaborate iterative method, where discretized systems are solved on finer and finer grids, in 
each step using information from the previous solution to choose only a subset of all possible 
constraints. This slows the growth of the number of constraints. (Details are described in 
Section 3.) With this step-wise mesh refinement scheme, as a proof of concept, we were able 
to obtain solutions on meshes with about 4500 points on each aperture using MATLAB and a 
standard LP solver. Using compiled computer language and more specialized solvers for optimal 
transportation problems, we expect that this number could still be substantially improved. 

There has been increased interest in numerical methods for optimal transportation problems 
in the last 10 years. Most work has concentrated on the Mongc- Ampere equation, which arises 
in the special case of optimal transport in K" with quadratic cost (also called optimal 
transport) [14, 2, 3, 7, 13, 23], although some authors have treated costs proportional to the 
distance (L^ optimal transport) [1]. These methods are based on fluid mechanical approaches 
[2] or varioiis finite difference approaches [3]. They are generally faster and allow for much 
larger mesh sizes than methods based on a discretization of the linear programming problem, 
but since they use the special structure of the Monge- Ampere equation, they are not directly 
applicable to more general cost functions or more general manifolds. 

One of the few papers that numerically addressed more general situations for optimal trans- 
port is [18]. A similar approach to the solution was taken there as we do in this paper, namely 
a discretization of the linear programming problem. (The authors in [18] used the primal for- 
mulation of the optimal transportation problem, whereas our approach is equivalent to using 
the dual formulation.) The authors noted that the number of variables is growing very fast 
with mesh sizes, making it impossible to solve the problem numerically even for moderate 
mesh sizes due to memory limitations. In fact, they noted that the software they used, the 
package Soplex [22], was designed to handle up to 2 million variables. This corresponds to 
mesh sizes of approximately 700 points on the input and output aperture in our problem if one 
employs a straightforward discretization scheme. The authors noted that for better results, one 
needs carefully designed programs. This is what we supply here for our problem: As explained 
above, our "mesh refinement" approach addresses this problem and makes it feasible to solve 
the problem for finer meshes. 

We note that our proposed algorithm does not make any a priori symmetry assumptions on 
the form the reflectors. It can also very easily be adapted for the numerical solution of other 
beam shaping problems for which a formulation using optimal transportation theory have been 
found, for example those in [9, 11, 20]. 

This article is organized as follows: In section 2, we recall the reformulation of the reflector 
construction problem as a linear programming transportation as given in [10] and fix some 
notation. We then describe a basic discretization scheme for the numerical solution and propose 
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Figure 1: Sketch of the reflector problem. The point source is located at the origin of the 
coordinate system. Note the conventions on coordinates as illustrated by the coordinate system 
in the lower left hand corner: The output beam propagates in the direction of the negative 
z— axis, and points in the plane perpendicular to the z— axis are denoted by the vector x G M^. 
Thus points in three dimensional space are denoted by (x, z). See Section 2 for more. (Adapted 
from [10].) 

an improved "mesh refinement" scheme in section 3. The following section 4 is devoted to 
numerical tests of this scheme. We first derive an explicit analytical solution to be used as 
synthetic data for our numerical work. Then we compute the solution numerically and analyze 
the error of approximation. Our scheme requires a sequence e^^-* , e*-^^ , ... of "constraint inclusion 
thresholds." We analyze different choices of this sequence. We conclude with a brief summary 
and propose future work. 

2 Formulation as Linear Programming Problem 

We now briefly review the notation for the problem posed in [15], as well as the result of [10], 
which reformulates the reflector design problem as a linear programming problem. 

Consider the configuration show in Figure 1. A point source at the origin O = (0,0,0) 
generates a spherical wave front over a given input aperture D contained in the unit sphere S^. 
The input beam has a given intensity distribution. By means of two reflectors, this wave front 
is to be transformed into a beam of parallel rays propagating in the direction of the negative 
z— axis. This output beam has to have a prescribed intensity distribution. The cross section of 
the output beam in a plane perpendicular to the direction of propagation is called the output 
aperture, and denoted by T. (Certain regularity conditions apply to D and T, see [10] for the 
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technical details.) 

Denote points in space by pairs (x, z), where x e is the position vector in a plane 
perpendicular to the direction of propagation and S R is the coordinate in the (negative) 
direction of propagation. See again Figure 1 for our convention on the direction of the z— axis. 
Points on the unit sphere will typically be denoted by m e 5^; their components are also 
written as m = (niajjm^) with jnij^p + ml = 1. 

We fix the output aperture in the plane z = —d. We will seek to represent the first reflector 
as the graph of its polar radius p(m) (for m G D), and the second reflector as the graph of a 
function z{x) (for z e T). (See Figure 1.) That is 

Reflector 1: Ri = {p(m) • m | m e -D}, Reflector 2: R2 = {(x, z{-k)) | x e T}. (1) 

The geometrical optics approximation is assumed. It follows from general principles of 
geometric optics that all rays will have equal length from (0,0,0) to the plane z = —d\ this 
length is called the optical path length and will be denoted by L. We define the reduced optical 
path length as £ = L — d. 

Oliker [15] showed that local energy conservation translates into a complicated partial 
differential equation of Monge- Ampere type for p(m). As noted in [15], the resulting equation 
is quite involved and a rigorous analysis of this equation seems very difficult. (See equation 
(59) in [15].) 

To amend this, the problem was reformulated in [10] as a linear programming problem, 
which makes a complete analysis possible, both concerning theoretical results on existence and 
uniqueness, and gives a method for practical computations. 

For this, the following function i4r(m, x), called the cost function in analogy with the theory 
of optimal transportation, plays an important role: 

£ — (m x) 1 - - 

'^("^'"^= 2^(^^-|x|^)a + m.) "4^ form=(m.,m.)eI?,xeT. 

In further preparation, the following two transformations are needed: 

Definition 2.1. Let z = z{x.) be a continuous function defined on T C R^. Then define the 
function 

^W-^-^^ forxef. (2) 

Definition 2.2. Let p = p(m) be a continuous function defined on D C with p > 0. Then 
define the function 

p(m) = -^^ + ^ , ^ ] — for m e D. (3) 

'^^ ' 2£ 2p(m) • (m^ + 1) ^ ' 

With these preparations, in [10], the following notion of a quasi-reflector pair and its as- 
sociated ray tracing map is used. (The term used in [10] is "reflector pair", but we use 
"quasi-reflector pair" here for clarity.) 

Definition 2.3. A pair (p, z) G C^q{D) x C{T) is called a quasi-reflector pair if p,z > and 



/3(m) = sup ( _j ^ ^(m, x) ) for m e D, (4) 
z{x)= sup ( —^K{m,x)] forxeT. (5) 
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Definition 2.4. Let (p, z) e C>o(-D) x C{T) be a quasi-reflector pair. Define its reflector map, 
or ray tracing map, as a set-valued map 7: 5 — >■ {subsets of T} via 




In [10], it is shown that 7(m) is in fact single- valued for almost all m e -D. Therefore 7 may 
be regarded as a transformation 7 : D ^ T. 

If {p,z) is a "quasi-rcflcctor pair" in the above sense, the following can be shown [10], 
justifying the choice of nomenclature: If physical copies of the corresponding surfaces (1) are 
made from a reflective material, then a ray emitted from the origin in the direction m € D 
will be reflected to a ray in the negative z— direction labeled by x = 7(m) S T. 

With this, the reflector problem can be formulated rigorously: 

Problem I. (Reflector Problem) For given input and output intensities /(m), m S D, and 
L(x),x e T, respectively, satisfying global energy conservation J^/(m)rf(T = Jj.L{x)dx, find 
a pair {p,z) G C>o(-D) x C{T) that satisfies the following conditions: 

(i) {p,z) is a quasi-reflector pair in the sense of Definition 2.3. 

(a) The ray tracing map 7: 5 ^ f satisfies 



Here da denotes the standard area element on the sphere S^. The second condition is local 
energy conservation. 

As we indicate below, the Reflector Problem can be reformulated in the following form: 
Problem II. Minimize the functional 



on the set Adm{D,T) = {(r,C) e C{D)xC{T) \ r(m) + C(x) > logiV'(m,x) for allinGD,xG 
f}. 



Problem II is the dual problem to the problem of finding a map P from D C to 
T C which minimizes the transportation functional P 1-^ K{ni, P{ni))da. Problem I 
and Problem II are equivalent, as expressed in the following theorem. 

Theorem 2.5. [10] Let {p,z) € C>o(-D) x C{T) be a quasi-reflector pair. Then (log p, log 0) G 
Adm(D, T). The following statements are equivalent: 

(i) {p,z) solves the Reflector Problem I. 

(a) (log /5, log 5) solves Problem II. 

Thus solving the reflector construction problem in Problem I is equivalent to solving the 
infinite dimensional LP problem in Problem II. Indeed, it can be shown that a solution as 
in Theorem 2.5 exists. (See Corollary 6.5 in [10].) In the remainder of this paper, we will 
concentrate on numerical solutions for this problem. 




for any Borel set t CT. 
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3 Numerical Schemes 



In the following, wc describe two scheme for solving Problem II numerically. First, let us list 
explicitly the given data on which the solution depends: 

(i) domains^ DC S'^,fCR'^ 

(ii) nonnegative integrable functions /(m), ui G D, and L(x), x G T, with Jq Ida = Jj, Ldx 

(iii) a reduced optical length (. 

(iv) the value p(mi) = pi for some fixed mi € D, where p is the radial function of the first 
reflector. ^ 

The flrst numerical scheme, Scheme 3.1 is a very straightforward discretization of Problem II 
via meshes on the domains D and T. Because of its simplicity and straightforwardness, we 
can informally call this a "brute force" technique. The problem with this scheme is that each 
pair of points from D and T gives rise to a constraint, and so the technique requires a lot of 
memory, so that only coarse meshes can be used. (This problem is very well known, see for 
example [18].) To address this issue, we propose a more sophisticated scheme. Scheme 3.2, 
which is an iterative scheme using finer and finer meshes, where in each iteration, information 
from the previous solution is utilized to reduce the number of constraints. 

The technique depends on choosing a parameter £, the constraint threshold, in each itera- 
tion. We then discuss the problem of choosing an appropriate sequence of such e. 

3.1 "Simple Discretization" Scheme 

The characterization of solutions by Theorem 2.5 immediately gives a practical solution scheme 
for solving the Reflector Problem I, namely by a straightforward discretization of the equivalent 
Problem II. (The same algorithm was used in [18].) 

Numerical Scheme 3.1. ("Simple Discretization") 

• Create a mesh in the input domain D by choosing sets rfi°\ d2°\ . . • , — ^^^^^ 

the interiors of any d^f^ and d^p are disjoint (for i ^ j), and l) = [J d^^ . (For instance, 

i=l 

rfi°\ d2°\ . . . , dj^j(o) may be a triangulation of D.) 

• Similarly, create a mesh in the output domain T by choosing sets , 4*^' , • • • , t^(o) ^ 

where the interiors of any t'f^ and ij*'^ are disjoint (for i ^ j), and T' = [J tf*^- 

• Choose sample points m^*^^ e d-"^ (for i = 1,...,M^^^) and x^*^' e tj^^ (for j = 

1, . . . , N^^'^). Here we may assume that m^^^ = mi is the same point as given in the 
data in (iv) above. 

^The rigorous result in [10] made the additional technical assumption (0,0, —1), (0,0, 1) ^ D. For practical 
purposes, these assumptions can often be dropped. 

^Without this constraint, the reflectors are not uniquely determined. This can be seen in the formulation of 
Problem II. Indeed, any transformation of the form ri-^r + c, C>-^z — c for constant c will leave the objective 
functional unchanged. 
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• Find the solution (^{rf^jfi^, {Cj°^}jii) of the following LP problem: 

Minimize ^ (^f) ^^ea (df^) + (f L (xf ) area (if) 

i=i j=i 

subject to r J*^' = log pi (where pi is as given in the data in (iv) above) , 
and rf ' + cj°^ > log K (mf , xf for » = 1, . . . , M^o) ; j = 1, . . . , iV^") . 

• Find the numbers ' ,i = l,..., M^o) and \ j = 1, . . . , Ar(o) such that rf ^ = log pf^ , 

and Cj*^^ = log Zj^^ ■ (This is straightforward by taking the inverse of the transformations 
given in Definitions 2.1 and 2.2.) 

Then pf'^ is an approximation for the true value p ^m^°^^ of the radial function of the 
first reflector, evaluated at the sample point m-'^' for i = 1,...,M^°). Similarly, Zj^^ is an 
approximation for the true value z (^f'^^ of the function describing the second reflector, eval- 
uated at the sample point x^'^' for j = 1, . . . ,N^'^\ (We arc using the superscript "(0)" here 
to distinguish this solution from additional iterative approximations we will obtain with the 
iterative scheme described below.) 

It is important to point out that this discretization also yields a discrctizcd version of the 
ray tracing map 7 (see Definition 2.4). Namely, for each index i = 1, . . . , there is at least 

one corresponding index j* , 1 < j* < N^^^ where the constraint corresponding to the pair 
is active, that is, where we have 

rf +Ci?=l0gi^(mf),x;.?). 

This means that the point nip is approximately the image of the point x^°^ under the ray 
tracing map 7. 

3.2 "Iterative Refinement" Scheme 

As noted in the introduction, one of the drawbacks of scheme 3.1 is that the constraint set in the 
linear programming problem in the penultimate step becomes large very fast with finer meshes, 
and the corresponding problem becomes too large to handle for standard LP solvers. We were 
not able to solve problems with more than very roughly 1000 mesh points on a standard PC 
with 4 GB RAM (M^o) « 500, N'^^^ « 500), which is arguably too small for many applications. 

We addressed this problem by developing an iterative scheme. First the problem is solved 
for a mesh with relativcy few sample points (say M^^^ = N^^^^ = 250). Then a finer mesh is 
chosen, and the previous solution is used to reduce the number of constraints. 

Specifically, we have the following scheme, depending on a number e^^^ > 0, which we call 
the "inclusion threshold," to be explained in detail below. 

Numerical Scheme 3.2. ("Iterative Refinement") 

• Use scheme 3.1 to find an initial solution {{rl^^}ff!i\{Cj^^}fli) for M^^^ sample points 
on D and N^°^ sample points on T, respectively. 
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• Create a mesh on D by choosing sets d^^^ ^2^^ . . . , cjj^^ij C D, where the interiors of any 
(fp and 

are disjoint (for i ^ j), and D = d\^\ Here M^^^ is chosen larger than 

i=l 

meaning that we have a finer mesh. 

• Similarly, create a mesh on T by choosing sets • • • : ^jvw — '^'^^'^^ interiors 

A,(i) 



of any t-^^ and tj^^ are disjoint (for i ^ j), and T = 



i=i 

• Choose sample points Txif"^ G df'^ (for i = 1,...,M'^^^) and x.^'' G 4"'^'' (for j = 

1, . . . , iV^^)). Here we may assume that m^^^ = mi is the same point as given in the 
data in (iv) above. 

• Interpolate the initial solution {rl^'^jfl^i^ to find an approximation for the first reflector. 
That is, by some interpolation method, find a continuous function r(m),m e D, with 

r(mf^) = rf^ for i = 1, . . . , M(°). 



• Similarly, interpolate the initial solution {Q }jLi to find an approximation for the 



.(0)-,jv(o) 

second reflector. That is, by some interpolation method, find a continuous function 
C(x),x e T, with 

C(xf)=cf fori = l,...,7VW. 

• Find the solution (^{rl^^}fii \ {Cj^hj^i) of the following LP problem: 

Minimize ^ ' / (m^) area (df + ^ L (xf area (t^) 
i=i j=i 

subject to r\^^ = log pi , 

and rf ^ + cj'^ > logiT (mf \xW) 

for those pairs (i,j) with r(mj^'') + C(Xj^'') ~ log if ^m^^\x^.^^^ < e^^\ 

• Find the numbers p\^\i = 1, . . . , M^^^ and z^^\j = 1, . . . , N^^'), such that rf'^ = logpf'\ 
and cj^' = log^j^^ 

The idea of the second scheme 3.2 is to solve the discretized LP on a coarse mesh first and 

then use this information to reduce the number of constraints needed for the LP on a finer 
mesh. To wit, the difference between schemes 3.1 and 3.2 is that in the discretized LP problem, 
in scheme 3.2 not all pairs of sample points from D and T, represented by pairs of indices {i, j), 
are included in the list of constraints. In fact, we would in principle only need to include those 
where the constraint is active, that is, where the corresponding inequality holds with equality. 
Of course, this information is not available a priori. Instead, we use the following heuristic: A 
constraint should only be included in the LP problem if it is "almost" active when we use an 
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interpolation of the solution on the coarse mesh as an approximate solution on the fine mesh. 
This is where the "inclusion threshold" e^^^ come into play. Namely, we only include those 
constraints for which the difference r(mj^^) + Ci^j^^) ~ ^ogK ^m^^\xj^^^ is less than 

the inclusion threshold £*^^^. (Here r(m) and C(x) are interpolations of the solution of the LP 
problem for the coarse mesh.) Clearly, increasing e^^^ means more constraints are included, but 
also potentially the solution of the LP represents a better approximation of the exact reflector 
pair. 

In our numerical tests, we found it advantageous to scale the functions /(m) and L{x) so 
that the approximations for the integrals Ida and Jj, bLdx yield exactly the same result. 

Note that the algorithm in scheme 3.2 can be iterated again. We can use the flrst iterative 
solution p[^\i = 1, . . . , M^^\ and Zj^\j = 1, . . . , N^^\ in Scheme 3.2 to obtain a second 
iterative solution p\^\i = 1, . . . , M^'^\ and Zj^\j = 1, . . . , N^'^\ and so on. 

3.3 Choice of thresholds for Iterative Refinement scheme 3.2 

For the iteration of scheme 3.2 described above, we have to choose a corresponding sequence 
of mesh sizes {(M("), iV(")), (M^i), A^(i))^ (Af(2), Ar(2))^ . _ (Here AfC^) denotes the number 
of sample points for the input aperture D, and A^f*^) denotes the number of sample points for 
the target aperture f , both in the k*^ iteration of the scheme.) We also need to choose a 
corresponding sequence of thresholds e*^^', e^^^ .... 

Let us assume that the sequence {{M^°'> , N^°^), {M^^\ N^^^), {M^^\ N^'^^), . . .} is given. 
This raises the question: How should the corresponding sequence of thresholds e^^\ s^^^ , e^^^ ,.. . 
be picked? 

This is a question of great practical importance. To wit, if the sequence {e*-'^^}/c is constant 
or decreases too slowly, then "many" constraints are included in each iteration, meaning that 
the LP problems will become large fast and we will run out of memory quickly after a few 
iterations. If the sequence {e^'°^}/c decreases too fast, then "few" constraints are included in 
each iteration. This could cause for instance that some index i (with 1 < i < M^''^) may not 
even be included in any of the constraints with any of the indices j (with 1 < j < A^^'^'), or 
vice versa. This would cause the LP problem to be unbounded and hence there would be no 
solution. (This behavior is illustrated in the practical tests summarized in Table 1; see also 
the discussion in section 4.) 

We can make a very rough estimate for a good choice of the sequence {e'^'^^jfe, assum- 
ing for simplicity that M^''^ w N^''\ Let us assume that each sample point m\''^ in the 
input aperture D is paired via the ray tracing map 7 with a unique point x^.^]^ w 7(m-''^) 
and vice versa. Thus there arc approximately M*^*^^ such pairs. Let r(m),m € D, and 
C(x),x S T, denote the function obtained by interpolating the solution of the (A; — 1)*'' it- 
erate, (^{r[''~^^}^^i ^' , {Cj''^^^}jili Then the set of all points in (m, x)— space that satisfy 
r(x) — z{x) — log K{ni, x) < e^*^' is approximately an ellipsoid for small e^'^h This can be seen 
by using the Taylor expansion around the pair (m^'^\ x^.^]^ ) in (m, x)— space. The volume of 

this ellipsoid is proportional to (e'-'^))^. Since each pair (mp'' , xj*^'' ) that satisfies this inequality 
corresponds to a constraint included in the LP of the fc*'* iteration, the total number of con- 
straints should be roughly proportional to M^'^^ ■ (e'*^^)^. So to keep the number of constraints 
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approximately constant throughout the iterations, based on these heuristics, one may choose 



£ 



where C is a constant. These heuristics are very rough and in practice the number of constraints 
is still increasing with M^'^^ albeit at a much slower rate than the case of constant e^^^ . In 
section 4.3, we describe an implementation of this idea. 

4 Numerical Tests 

To test the validity of the numerical scheme; described above, we used it on a case where the 
solution is known in analytic form. In the next subsection, we first describe this special analytic 
solution, and then discuss our results. 

4.1 An Analytic Solution 

In order to obtain a data set where an explicit analytic solution is known, we consider a special 
given configuration of reflectors, and then solve the "forward" problem of determining the 
output intensity this system produces for a given input intensity. 

4.1.1 Construction of a pair of reflectors i?2 

Consider the following set of two reflectors R\ and i?2, sketched in Figure 2: Let a = (a, 6, c) 
be a given point in R^, and let i? > and a > be two positive numbers with i? > |a| = 
Vo^ + 6^ + . Now let the first reflector i?i be the boundary of a spheroid with foci at the 
origin O and at a and with major diameter i?/2. (Thus for each point P on the sum of 
the distances from each of the two foci PO + Pa equals i?.) Let the second reflector i?2 be the 
boundary of a paraboloid whose main axis is the negative z— axis and whose focus is at a with 
focal parameter 2q. (Thus for each point p = {x,y,z) on R2, the sum of the distance to the 
focus and the shifted 2;— component Pa + {z — c) equals 2a.) 

Note that by definition of the two reflectors, any cone of light rays emitted from the origin 
O will be transformed into a beam of parallel rays traveling in the direction of the negative 
2:— axis. Indeed, a ray emitted from O will be reflected off Ri towards the focus a. This ray 
will then be reflected in the direction of the negative z— axis by R2. See again Figure 2. 

It is not hard to find explicit expressions for the two reflectors. If we write i?i = {p(ni) • 
m I m e 5^} and R2 = {(x, z{x)) | x = (x, y) G M^}, then we have the following expressions: 

R2 _ |a|2 

radial function p for Ri : p(m) = — r (6) 

and 

equation for R2 : — 4:aq = |p|^ — 4a^ (7) 
where we used the shifted coordinates (p, q) = {x, y, z) — (a, b, c) = (x, z) — a. 
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X = /(m) 



Reflector /?2 
(paraboloid) 



(x, z(x)) 



p(m)m 




a = {a,b, c) 



Reflector /?j 
(spheroid) 



^ z 



Figure 2: Sketch of the reflectors Ri and i?2- is the boundary of an eUipsoid whose foci are 
the points O and a. R2 is the boundary of a paraboloid whose focus is a. The sketch shows 
a two-dimensional cross section. A ray given by the direction m e 5^ will be reflected by i?i 
and then by R2 to a ray in the direction of the negative 2— axis. 
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4.1.2 Ray Tracing Map 7 

One can now determine the ray tracing map 7 corresponding to the reflector pair Ri, R2. See 
again Figure 2. 

Let m e S*^ be a given direction. Thus m = (nix, mz) with nix G and |mxp + m1 = 1. 
To determine the image 7(m), consider the ray emitted from the origin in the direction m. The 
ray will encoimter the first reflector at the point p(m) • m and will then be reflected towards 
the second focus a. Thus the reflected ray can be parameterized by a + A • (a — /9(m) • m), 
where A is a parameter. With equation (7), this yields that the point where the reflected ray 
hits the second reflector is a + A* • (a — p • m) where 

= {R-pr-l-pm.r " ^'-'""'^^ = R + c-pil + m.y 

The projection of this point to the xy— plane is the value of the ray tracing map. Thus we have 
the following explicit formula for the ray tracing map: 

^^"^^ = ( 6 ) + i? + c-p(m)(l + m.) (( 6 ) , m= (mx,m.) e S^. (8) 

4.1.3 Solution of the Forward Problem 

We can now solve the forward problem: Given the reflector pair i?i,i?2 as deflned above, an 
input aperture fj C 5^ and an intensity distribution /(m), m S f2, flnd the output aperture 
T C M2 a,nd the output intensity L(x), x € T. _ 

The output aperture is simply the image of fl under the ray tracing map 7: T = To 
flnd the induces intensity L{x), use the deflning property 

/ I{m)da^ I L(x)dx (9) 

J J 7(0;) 

for all Borel sets a; C O. This is an energy balance equation. 

The above integral equation allows us to find an explicit expression for i(x) given J(m), 
or vice versa. We consider for simplicity the case that Q. is contained in the left hemisphere 

= 5^ n {(a;, y,z) & B?\z < 0}. We can then use coordinates 

t: {{m^.my) €R^\ml + ml <l} ^ S'i, (10) 
{x,y) ^ (m^,my,-^l- ml - mlj . (11) 

In these coordinates, the standard measure on 5^ is given by da = ^"^^ ^ ^here 



iriz = — -i/l — ml — my. Using this, (9) is then equivalent to the equation 



/(mx, m,) = i(7(T(mx))) J(7 o r). (12) 



HereJ(7or) = |det( ^ ^) 



denotes the Jacobian of the map 7 or, With the help 
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of a computer algebra system like Mathematica, this can be evaluated explicitly as 

J(7 o r) = '- 2 ' (13) 

-m, (2{c + R)in^-(^ ° ^ - {1 + m,)\{a,b)\^ - {c + R)^{1 - m,)^ 

where again ruz = — y^l — — rriy. 
4.2 A data set 

We can now construct a particular data set using the results of the previous sections. We pick 
a = b = 0, and c = —0.4, a = 1, R = 1.3, and consider the input aperture 

^ = {{rrix, my, m^) C S"^ \ ^ml + ml < 0.8, < 0}, 

which is a spherical cap centered at the point (0, 0. —1) . Using the results from the previous 
subsections, we can now find the output aperture in a straightforward manner: 

f = {{x, y)eM.^\ ^/x^ + y2 < 1.8888888889}. 

We choose a constant output intensity 

L(x) = 1 for X e T. 

Using the relation (12), this gives the input intensity 

, 14.2716049383 , , ^ ^ 

i(m) = 2 — for m = (mx,my,m2) e S2. 

(1 - Tflz) 

The two reflectors are given as follows: 

0.765 

radial radius p for R\ : p(m) = — — — — (14) 

^ 1.3 + 0.4m2 ^ ^ 

and 

equation for R2: z = -0.25 (x^ + y"^) + 0.6. (15) 
The corresponding reduced optical path length is 

^ = i? + 2a = 2.9. 

This follows from the construction of the two reflectors and the properties of ellipsoids and 
paraboloids. The two reflectors are shown in Figure 3. 



^The data do not represent a physically possible set of reflectors since there would be blockage. However, 
the example is useful as a proof of concept, and the numerical scheme itself is completely independent of the 
shape of the apertures or any a priori symmetries of the problem. 
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Figure 3: Plots of the reflectors given by (14) (left) and (15) (right) in the explicit data set 
given in section 4.2. 



4.3 Results 

Using the data set from subsection 4.2, we conducted a series of numerical tests to establish the 
validity of the numerical scheme 3.2. Since we have explicit analytical expressions for the two 
solutions p{m) and z{x) describing the two reflectors as given in (14) and (15), respectively, 
we were able to compute the error of approximation directly. (See also Figure 3 for plots of 
these reflectors.) 

The iterative scheme 3.2 was implemented in MATLAB with the package distmesh [17] 
for mesh generation in each iteration step and lp_solve [4] for solving the resulting sparse LP. 
Computations were performed on an Intel 15 Dual Core/2. 53GHz with 4GB RAM. Our primary 
aim is to demonstrate the practicability of the proposed scheme. While there are several ways 
in which the implementation could be made more efficient, for instance by using a compiled 
computer language or a more specialized LP solver for transportation problems, we believe 
that our results show that the proposed scheme is easy to implement and makes it possible to 
use much finer meshes than a simple discretization scheme. 

The mesh generation algorithm employed by the software we used is based on a relaxation 
scheme of forces in a truss structure [17]. For this, user input is the desired edge length /i, 
not the number of mesh points M (for the domain Cl) or N (for the domain T), respectively. 
However, it is not hard to see that the relation between the average edge length h and the 
number of mesh points M on D obeys 




Indeed, the total area A is proportional to the number of mesh points M times the area of one 
triangle of a Delaunay triangulation. The area of a triangle is proportional to the square of 
the average edge length h. The same relation holds of course for the number of mesh points N 
on T. (In Tables 1 and 2 in this section, we show the number of mesh points M, N instead of 
the edge length h, as we believe the former are more informative.) 

We chose a sequence of desired edge lengths h^'^\ h^^\ . . .. Then, based on the heuristics 
given above and in section 3.3, we determined the corresponding constraint thresholds with 
the formula 
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eW=C. (fc = 1,2,3,...), (17) 

where C and a are constants. We tested several values for the proportionality constant C and 
the exponent a. (The heuristic arguments in section 3.3 yields a = 1, but we also tested other 
values.) Some results arc summarized in Tabic 1. As indicated in section 3.3, "large" values for 
C mean that "many" constraints are included in the LP in each step, which potentially means 
good approximations, but also high memory usage, and thus we may run out of memory after 
a few iterations. In contrast, "small" values for C improve memory usage, but may mean that 
the numerical results may be less accurate approximations of the exact results, or indeed the 
problem may become unbounded, as typically seems to happen in practice. (See Table 1.) A 
similar logic holds for the exponent a: For larger a, the threshold e*^*^^ decreases faster, leading 
to less memory usage, but the danger of the problem becoming unbounded. 

Table 1 confirms and illustrates these principles. At each iteration, the LP is either un- 
bounded or there is a solution. Interestingly, if there is a solution, it seems largely independent 
of the number of constraints, as indicated by the fact that the error of approximation varies 
very little. (See fourth column in Table 1.) This makes some intuitive sense as the most im- 
portant issue in each iteration is really that we have included the active constraints. As long 
as all active constraints are included, we get the same solution. If however we did not include 
any of the constraints for a certain mesh point, the problem becomes unbounded. 

In fact. Table 1 illustrates the following feature of the algorithm, which is intiiitively quite 
obvious, although we do not give a formal proof: At each iteration step, there is a critical value 
£crit such that if e < Ecrit, then the LP is unbounded because no constraints involving a certain 
point are included. If £ > Ecrit, the LP has a solution. As indicated above, it also appears that 
if there is a solution, it is likely typically close to the solution of the "full" problem (that is, 
the problem including all possible constraints). This leads to a possible alternative algorithm 
for choosing the thresholds e*^^',e(^\ . . . : Instead of choosing each c^*^-* with a pre-determined 
formula, use the corresponding critical value Ecrit in each iteration step. This critical value 
can be found by simple trial and error. While this idea is promising, it is possible that such 
a scheme may lead to a much slower decay of the maximum error, or even an increase of the 
error. We did not test this idea in this article, but it would be interesting to analyze it further. 

To find out the maximum mesh size on which a solution could be obtained without running 
out of memory or arriving at an unbounded problem, we "calibrated" the values of C and a. 
The process is summarized in Table 2. 

We obtained best results (that is, largest mesh sizes) for the run with C = 1.7 and a = 1. 
Table 3 now summarizes the results for this specific run in more detail. The table shows the 
errors between the numerical computed approximations and the exact formulas as given in (14) 
and (15). Figure 4 shows plots of the approximation errors for the reflectors as functions on 
the input and the output apertures. 

These results are a good indication that the scheme converges, and that the error decreases 
proportional to M~", or equivalently proportional to /i^", where roughly 0.5 < a < 1. (The 
results from Figure 1 appear to indicate a ~ 1, but in other data, we also encountered results 
where a seemed to be closer to 0.5.) The results also illustrate the practicability of the scheme. 
Note that the number of mesh points was increased by a factor of about 16 from iteration 1 
to iteration 6. In fact the problem in iteration 6 would be impossible to solve with a simple 
discretization scheme due to the size of the problem. In iteration 6, we were able to reduce 
the number of constraints used to only about 15% of the number of constraints used for the 
"naive" simple discretization scheme 3.1. 
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c 


a 


sequence of mesh points 


notes 


max no. constr. 


% of possible constr. 


1 


1.1 


284;455;724 


unbounded on 724 


N/A 


N/A 


1.2 


1.1 


284;455;724;1147 


unbounded on 1147 


N/A 


N/A 


1.3 


1.1 


284;455;724;1147 


Max err 0.0018575 


222220 


16.9% 


1.5 


1.1 


284;455;724;1147 


Max err 0.0018521 


258185 


19.6% 


2 


1.1 


284;455;724;1147 


Max err 0.0018521 


342846 


26.1% 


5 


1.1 


284;455;724;1147 


Max err 0.0018521 


713571 


54.0% 


0.8 


1 


284;455;724 


unbounded on 724 


N/A 


N/A 


0.9 


1 


284;455;724;1147 


unbounded on 1147 


N/A 


N/A 


1 


1 


284;455;724;1147 


Max err 0.0018521 


226231 


17.1% 


1.5 


1 


284;455;724;1147 


Max err 0.0018521 


340038 


25.7% 


1 


1.5 


284;455 


unbounded on 455 


N/A 


N/A 


2 


1.5 


284;455 


unbounded on 455 


N/A 


N/A 


3 


1.5 


284;455;724;1147 


unbounded on 1147 


N/A 


N/A 


3.5 


1.5 


284;455;724;1147 


unbounded on 1147 


N/A 


N/A 


4 


1.5 


284;455;724;1147 


Max err 0.0018521 


224140 


16.9% 


5 


1.5 


284;455;724;1147 


Max err 0.0018521 


282415 


21.4% 



Table 1: Table summarizing the resuhs for rmis of the iterative scheme with different values for 
the constants C and a given in (17). The third column shows the sequence of mesh points M*^*^) 
for D. (So for instance, in the first row, in the initial discretization, there were 284 points. In 
the first iteration, there were 455 mesh points. The second iteration had 724 points, but the 
problem was unbounded.) The third column notes if the problem became unbounded before 
reaching 1147 mesh points. If the total run was completed up to 1147 mesh points, the third 
column gives the maximum norm Hpnum — Pana||oo of the error between the resulting numerical 
solution Pnum and the analytic solution Pana- The fifth column shows the number of constraints 
for the case of 1147 mesh points, if the run reached this maximum number. The last column 
indicates the corresponding percentage of number of constraints relative to the total number 
of constraints in the "brute force" scheme 3.1. (So for example in the last row, there were only 
21.4% of all possible constraints used in formulating the last LP, thus eflfectively reducing the 
number of constraints needed by 78.4%.) 
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c 


Oj 


no itPTatioTm 


no TTIP9V1 ■nnint^ 




no poTisf.r 


% of ■nos9iV)lp POTi^f.r 


n R 


\ 


2 


724 


Mfl-sr prr 00148 


101174 


X C/.'rt /(J 


n R 


\ 




1148 


1 1 Tl r»/~n 1 Tl H aH 


1 78976 


1 3 6% 

XO.U /(J 


1 


1 


3 


1148 


Maif prr 001^,^21 

XVJ.CUA. C'J. J. Vy • V-' Vy J- O ^ ^ J- 


226231 


17 2% 


1.2 


1.1 


3 


1148 


Mfl-sr orr flOl 8575 


203535 


15.5% 


-L .0 


1.1 


3 


1148 


May prr 001857^ 


922220 


16 9% 




1.1 


3 


1148 


Mav prr fl 0018591 


258185 


1 9 6% 


2 


1.1 


3 


1148 


Max err 0.0018521 


342846 


26.1% 


1 


1 


4 


1824 


Unbounded 


446239 


13.5% 


1.5 


1 


4 


1824 


Max err 0.00131 


688163 


20.7% 


1.6 


1 


4 


1824 


Adax err 0.00131 


733894 


22.05% 


1.5 


1 


5 


2882 


Max err 0.00069 


1373976 


16.5% 


1.6 


1 


5 


2882 


Max err 0.00069 


1472275 


17.7% 


1.5 


1 


6 


4536 


Unbounded 


2650000 


12.8% 


1.6 


1 


6 


4536 


Max err 0.00067 


2856627 


13.88% 


1.7 


1 


6 


4536 


Max err 0.00067 


3057070 


14.86% 


1.6 


1 


7 


7130 


Out of memory 


N/A 


N/A 


1.7 


1 


7 


7130 


Out of memory 


N/A 


N/A 



Table 2: Table illustrating the process of calibrating C and a from (17) for maximum mesh 
sizes. Each row represents a different run of the iterative numerical scheme. The third column 
is the number k of iterations in the run and the fourth is the number M^^^ of points in the 
mesh for D in the last iteration. The fifth column notes the result of the run - if a solution 
was found, the maximum error is given. If not, the linear program either became unbounded, 
or the computer ran out of memory. The last two columns are the number of constraints used 
for the last iteration and the percentage of total possible constraints (compared to the "brute 
force" scheme 3.1). 



k 


h 


M 


N 


Constraints 


Max error Ri 


L2 error Ri 


Max error i?2 


L2 error R2 


1 


0.12 


284 


278 


78,952 


0.0048 


0.00143 


0.008 


0.0021 


2 


0.096 


455 


450 


85942 


0.0022 


0.00076 


0.0047 


0.0014 


3 


0.0768 


724 


721 


183308 


0.00148 


0.00056 


0.0039 


0.0012 


4 


().()()144 


1148 


114() 


381971 


().()()12 


().()()()39 


o.ooLsr) 


().()()()44 


5 


0.049152 


1824 


1810 


779053 


0.00060 


0.00021 


0.0013 


0.00033 


6 


0.0393216 


2882 


2879 


1,568,533 


0.00059 


0.00019 


0.00069 


0.00016 


7 


0.0314573 


4536 


4525 


3,057,070 


0.00045 


0.00010 


0.00067 


0.00027 



Table 3: Error for a run of the iterative scheme 3.2 with C = 1.7 and a = 1 in (17). First 
column: iteration number. Second solumn: average edge length h. Third and fourth columns: 

number of mesh points for the domains of reflector 1 and reflector 2, respectively. Fifth column: 
resulting number of constraints in the LP. The scheme stopped after iteration 7 because there 
was not enough memory available for the next iteration. The maximum error as a function 
of total mesh points A'tot = N + M decayed in the form Smax ^ (A'tot)" with a w —0.82 for 
reflector 1 and a ~ —0.95 for reflector 2. (These values for a were obtained through least 
squares curve fitting.) 
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Figure 4: Error surface plots for iteration 7 in Table 3 for reflector 1 (left) and reflector 2 
(right). For better visibility, the plots were obtained by interpolating the error at each mesh 
point to a continuous function. Note that the error distribution appears to be relatively 
uniform, although the "spikes" are noticably arranged in triangular patterns and in a fashion 
somewhat reminiscent of spokes on a wheel. This may be an artifact of an iterated error in 
conjuncture with triangular meshes. 

5 Conclusion and Future Work 

We proposed two numerical schemes for solving an infinite dimensional optimal transportation 
problem arising in reflector design: A straightforward discretization and an improved iterative 
scheme, which uses knowledge of the previous solution in each step to reduce the number of 
constraints. This scheme is easily adapted to similar transportation problems arising in beam 
shaping problems, e.g. in [9] and [11]. We showed that this new scheme is easy to implement 
and makes it possible to solve the problem on much finer meshes. 

There are a number of possible directions for future research. We did not rigorously prove 
that the scheme converges, although we strongly expect that it does. In fact the decreasing 
error shown in Table 3 gives evidence for this, at least in our example. It would be valuable to 
have a rigorous proof of convergence. 

Another direction for further investigation is to change the selection algorithm for the 
threshold values e^^\ e^^^ ... as suggested in section 4.3: It is a straightforward conjecture that 
in each step, there is a critical threshold inclusion number such that the problem is unbounded 
for values of e below that number. One way to select a value of e is to pick a value for e*^*^) 
close to the critical value £crit. It would be interesting to flesh out this ideas and analyze the 
results more closely. 

Another avenue of research is to see whether our proposed solution algorithm can be made 
parallizable. 

It would also be worthwhile to investigate whether the existing PDE-based schemes for 
solving the transportation cost with quadratic costs [2, 3, 7, 23] can be adapted to the problem 
at hand. 
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